clc; clear;

%% 参数设置
alpha_list = linspace(2, 3, 20);    % alpha参数范围
q10_list = linspace(0, 1, 20);      % q10参数范围

A_template = [-1, 0, 0;
               0, 0, 0;
               0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

dt = 0.05;            % RK4步长
T_total = 600;        % 总积分时间
T_trans = 500;        % 暂态时间
N_steps = round(T_total/dt);
N_trans = round(T_trans/dt);

delta0 = 1e-8;        % 微扰大小

LE_vec = zeros(length(q10_list)*length(alpha_list),1);

fprintf('开始计算最大李雅普诺夫指数 (RK4积分)...\n');
parfor idx = 1:length(alpha_list)*length(q10_list)
    [i_q10, i_alpha] = ind2sub([length(q10_list), length(alpha_list)], idx);
    alpha = alpha_list(i_alpha);
    q10 = q10_list(i_q10);

    A = A_template;
    A(1,1) = -alpha;

    x0 = [1e-6; 0; 0; q10];

    LE = estimate_max_LE_RK4(@(x) mCNN_single(0,x,A,B), x0, dt, T_total, T_trans, delta0, N_steps, N_trans);

    LE_vec(idx) = LE;
end
fprintf('计算完成。\n');

LE_map = reshape(LE_vec, [length(q10_list), length(alpha_list)]);

%% 绘制热力图
figure;
h = imagesc(alpha_list, q10_list, LE_map);
set(gca,'YDir','normal');

colorbar;

% 设置NaN透明
alpha_data = ~isnan(LE_map);
set(h, 'AlphaData', alpha_data);

xlabel('\alpha');
ylabel('q_{10}');
title('最大李雅普诺夫指数');
grid on;


%% --- Wolf轨道扰动法计算最大LE函数（RK4积分版） ---
function LE_max = estimate_max_LE_RK4(f, x0, dt, t_total, t_trans, delta0, N_steps, N_trans)
    x1 = x0;
    x2 = x0 + delta0 * [1; 0; 0; 0];  % 固定扰动方向
    
    sum_ln = 0;
    count = 0;
    
    for i = 1:N_steps
        % RK4积分主轨道和扰动轨道
        x1 = RK4(f, x1, dt);
        x2 = RK4(f, x2, dt);
        
        dist = norm(x2 - x1);
        if dist == 0
            dist = 1e-16;
        end
        
        if i > N_trans
            sum_ln = sum_ln + log(dist / delta0);
            count = count + 1;
        end
        
        % 归一化扰动距离
        diff_vec = (x2 - x1) / dist;
        x2 = x1 + delta0 * diff_vec;
    end
    
    LE_max = sum_ln / (count * dt);
    
    % 绝对值大于10的设为NaN
    if abs(LE_max) > 10
        LE_max = NaN;
    end
end
